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Abstract The diffusion of finite-size hard-core interacting particles in two- or 
three-dimensional confined domains is considered in the limit that the confine- 
ment dimensions become comparable to the particle's dimensions. The result is 
a nonlinear diffusion equation for the one-particle probability density function, 
with an overall collective diffusion that depends on both the excluded-volume and 
the narrow confinement. By including both these effects the equation is able to 
interpolate between severe confinement (for example, single-file diffusion) and un- 
confined diffusion. Numerical solutions of both the effective nonlinear diffusion 
equation and the stochastic particle system are presented and compared. As an 
application, the case of diffusion under a ratchet potential is considered, and the 
change in transport properties due to excluded-volume and confinement effects is 
examined. 

Keywords Brownian motion ■ Fokker-Planck equation • Diffusion in confined 
geometries • Entropic effects • Stochastic simulations 



1 Introduction 

Transport of material under confined conditions occurs throughout nature and 
applications in industry. Examples include the transport of particl es in biolog ical 
cells, such as ion channels that conduct ions acr oss the cell surface (iHille 2001) or 



intracellular cargo alo ng microtubule filaments (jAlberts et a l l2002HKlumpp et all 



2005h . and in zeolites (|Keil et alll2000l) . Similarly, confinement can be i mportant in 



the d iffusion of cells themselves (e.g. blood cells through microvessels. lPries et all 
Il996l ) and surface diffusion on the cell me mbrane, which is usually crowded with 



fixed and mobile obstacles (jNicolau Jr. et al .2007). Moreover, recent advances in 



nanotechnol ogy have allowed the developme nt of synthetic nanopores and microfiu- 
idic devices (jHanggi and Marchesonill2009l) . which can be used for the sensing of 
single particles (such as small molecules, organic polymers, proteins, or enzymes) 
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and for studyi ng chemical reactions, biomolecular rec ognition, and interactions at 
the nanoscale (|Dekkedl2007l :[ Howorka and Siw A common feature in these 



applications is the interplay between the particle motion (usually noisy) and the 
geometric constraints. An additional factor comes into play if the system contains 
a collective of interacting particles rather than an individual particle. The way in 
which these characteristics combine to produce the global behaviour is a crucial 
factor in the understanding of such systems. In many respects, these transport 
phenomena can be studied in t erms of the canoni cal problem of geometrically 
constrained Brownian dynamics (jBurada et a 



When considering a theoretical model of particle diffusion in confined envi- 
ronments, there are three important modelling decisions to make. First, one must 
decide on the most appropriate representation of the particle diffusion and inter- 
actions (with other particles and the confining walls) . For exampl e, a common ap- 
proac h is to use a lattice-based random walk model with exclusion ( Plank and SimpsonI 



l2012l) . that is, to assume that the motion of particles is restricted to taking place 
on a lattice and that any attempted move to an occupied site is aborted. An alter- 
native approach is to consider a lattice-free random walk, in which the individual 
particle movements are not restricted to a lattice. It this case, excluded-volume in- 
teractions can be taken into account by assuming particles are hard spheres which 
cannot overlap each other, t hus considering a Brownian motion of hard spheres 
( Bruna and Chapman|[2012bl) . While in some cases a lattice-based model is more 



suitable f or the part icular application, in general the lattice-free approach is more 



realistic ( Plank and Simpson ,20l3) and the choice of an on-lattice model is for 



technical convenience only. 

The second modelling decision concerns the level of description, that is, whether 
to use an individual-based model or a population-based model. In the first case, the 
system of diffusing and interacting particles is represented with a stochastic model 
that describe the dynamics and interactions of each particle explicitly. This is typ- 
ically a computationally intensive approach, involving many statistically identical 
realisations of the stochastic simulation to develop insight into the population-level 
dynamics. In contrast, the population-based model consists of a continuum descrip- 
tion of the system in the form of a partial differential equation (PDE) for the popu- 
lation density of individuals. The continuum model tends to be easier to solve and 
analyse and can be particularly useful when, for large systems of interacting parti- 
cles, discrete models become computationally intractable. However, the challenge 
is to predict the correct PDE description of a given system of interacting particles, 
and, as a result, many population-based models are described phenomenologically 
at the continuum level rather than derived from the underlying particle transport 
process. For example, while it is well-understood that a non-interacting Brown- 
ian motion is associated a linear diffusion PDE at the population-level, it is not 
so straightforward to predict how excluded-volume interactions at the d iscret e 
level emerge in the PDE model. As pointed out in IPlank and Simpso 3 (|2012l) . 



the ability to represent mathematically both the individual-level details and the 
population-level description of a stochastic particle system is important because 
many experimental observations involve data at both levels for the same system. 
As a result, if we are to use both the individual-based and the population-based 
models of the same system, the link between the two must be fully understood to 
ensure that both models are consistent with each other. 
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Finally, the third consideration has to do with the way confinement is included 
in the model. It is important to note that the idea of confinement is inevitably 
relative the particle's characteristic size. A common approximation of confine- 
ment situations which is applicable when the particles are much smaller than the 
channel width is to ignore steric interactions between particles (assuming they are 
simply points) and only consider the geometric effects of the confining environment 
( Burada et a 1200^. For example, the diffusion of point particles in a (narrow) tube 
of varying cross-section can be approximated by a n effective one-dimensional diffu- 
sion e quation known as the Fick-Jacobs equation (|jacobslll967l :lR eeuera and Rubil 
l200lh . Another example in which particle interactions are omitted can be found 
in the Brownian ratc het models of molecular motors ([Munoz-Gutierrez et al 2 0121 : 
lEichhorn eraill2002[ ). which take the form of a one-dimensional diffusion under a 
periodic potential and tilting force. 

The opposite limit is single-file diffusion ( Henle et alll2008h . in which the finite- 
size of particles is taken into account but the confinement is so extreme that parti- 
cles cannot diffuse past each other (imagine a channel of width equal to the diame- 
ter of particles). Mathematically this problem is modelled as a one-dimensional do- 
main with hard-c ore interacting particles (hard rods) and has been wi dely studied 
(see, for example. iLizana and Ambiornssonll2009HBodnar and Velazquez 2005h . 

Both of these limits are extremes. The distinguished limit in which the finite- 
size interactions are important but the confinement is not so extreme that particles 
cannot pass one another has rec eived little attentio n; one notable exception is the 
exclusion process on a lattice in lHenle et all (200^. 



1.1 Aim of this paper 

This work introduces a theoretical framework for studying particle diffusion pro- 
cesses in confined environments. Rather than attempting to answer a particular 
question related to one of the applications presented earlier, here we are interested 
in developing a technique to tackle the common first steps in any of these such 
problems. Following the three considerations outlined above, we are interested in 
a lattice-free approach, in deriving the population-level model systematically from 
the individual-based model, and in an intermediate level of confinement. To this 
end, we consider the evolution of a system of TV identical hard spheres in a confined 
domain, in the limit that the confinement dimensions become comparable to the 
particle dimensions. In this setting, the finite size of particles is important not only 
for particle-particle interactions, but also for interactions with the domain walls. 
We consider in particular three confinement scenarios: a two-dimensional channel, 
a three-dimensional square channel, and two close parallel plates. However, since 
our approach is systematic, our model can be extended to other geometries. 

The key idea is that the system will reach equilibrium in the confined directions 
quickly, leading to an effective diffusion of reduced dimension in the unconfined 
directions only. With this in mind, the solution procedure consists of two steps: 
first, to reduce the model of A'' interacting particl es to a model for the evolutio n 
of the one-particle marginal density, as we did in iBruna and Chaomanl (|2012bl) : 



and second, to reduce the resulting model from a d-dimensional confined domain 
to an effective one-dimensional axial model in the case of a narrow channel, or to 
an effective two-dimensional planar model in the case of parallel plates. 
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1.2 Plan of this paper 

The work is organised as follows: in the next section we will introduce the problem 
setup, illustrate how the problem simplifies in the case of point particles and 
present the main result of this work, a population- level PDE model for the diffusion 
of hard spheres as a function of a confinement parameter, given by equation (|10p . 
In the third section we examine how our model interpolates between the different 
limiting cases of confinement. In Section 4 we explore numerical solutions of our 
PDE model and compare them with stochastic simulations of the particle-based 
model and numerical solutions of the limiting models. Finally, the fifth section will 
be devoted to the derivation of (jlOp for a two-dimensional channel. 

2 The model 

2.1 The setup: drift-diffusion in confined geometries 

We consider a population of N identical particles diffusing in a bounded domain 
n C M*^ (d = 2,3), interacting with each other and the domain walls with a 
repulsive hard-core potential, and in the presence of an external force. We work in 
the dimensionless problem by scaling space with a typical unconfined dimension 
L, time with /Do where Dq is the constant molecular diffusion coefficient, and 
force with jDo/L where 7 is the frictional drag coefficient. We assume particles 
are spherical with nondimensional diameter e ^ 1. 

Assuming the overdamped limit, the stochastic dynamics of the system is de- 
scribed by a set of stochastic Langevin equations 

X,(f) =f(X,(t))t + %/2W,(t), i = l,...N, (1) 

where X.i{t) G f} denotes the centre of particle i at time t > 0, f is the dimension- 
less external force (or drift) and W.^ are N independent d-dimensional standard 
Brownian motions. We note that by writing f(Xi(t)) we are assuming that the 
force acting on the ith particle only depends on its own position, thus excluding 
forces sucli as the electromagnetic force which would depend on tlie positions of 
all the particles X = (Xi, . . . ,Xjv). We suppose that the initial positions Xi(0) 
are random and identically distributed. Note that, because of the finite size of 
particles, we have tire set of constraints ||Xi — Xj || > e for i ^ j, so that the system 
of SDKs (P) is coupled. 

The Langevin system ([1]) is equivalent to the Fokker-Planck equation for the 
joint probability density P{x, t) of the A'' particles to be found at the position 
X = (xi, . . . ,xjv) e at time t, given by 

^ (f , t) = Vs- [VgP - F{x)pj , (2a) 

where Vg and Vg ■ respectively stand for the gradient and divergence operators 
with respect to the A^-particle position vector x and F{X) = (f (5), . . . , f (a;)) is the 
total drift vector. Because of excluded-volume effects, the domain of definition of 
^ (or configuration space) is not but its hollow form = \ Be, where 
Be = {x £ : 3i ^ j such that ||xi — Xj || < e} is the set of all illegal configurations 
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(with at least one overlap). On the contact surfaces dfi^ we have the reflecting 
boundary condition 

0= [VsP-F(f)P] -n, (2b) 

where n e gd.N-i jjgj^Q^gg ^^jjg ^jjj^ outward normal. Finally, since the particles are 
initially identically distributed, the initial probability density P{x,0) = Po{S) is 
invariant to permutations of the particle labels. The form of Q then means that 
P itself is invariant to permutations of the particle labels for all times. 

We suppose that J7 is a confined domain, with k < d confinement dimensions 
which are comparable to e. We introduce de = d — k as the effective dimensionality 
of the problem. In particular, we shall consider the following cases: 

o (NC2) Two-dimensional narrow channel (d = 2, fc = 1 and de = 1): 

^= [-il] X hf.f]- (3a) 
o (NC3) Three-dimensional narrow channel (d = 3, fc = 2 and de = 1): 

r2=[-i,i]x[-f,f]x[-f,f]. (3b) 

o (PP) Two parallel plates (d = 3, fc = 1 and de = 2): 

r?=[-i,i]x[-i,i]x[-f,f], (3c) 

where H = 0{e) is the confinement parameter. We note that H > 0, with H = 
allowed since J7 is the volume available to the particles' centres. In the case of a 
narrow-channel, when H < I particles cannot pass each other. We assume that 
the volume fraction is small; since \n\ = 0(e'=) this implies that Ne''-" < 1. 



N particles, d dimensions 


Ti 


1 particle, d dimensions 


P(xi, . . . ,xjv,t) 




p(xi,t) 










1 particle, 1 dimension 




Pe(£l, i) 



Figure 1 Schematic of the problem solution steps for de = 1 [narrow-channel cases (NC2) 
and {NC3)]. The goal is transformation T, to obtain an effective one-dimensional equation 
along the channel for the marginal density of one particle. We achieve this with the combined 
steps 71 followed by 72- 



The high-dimensional diffusion problem ([2]) will be reduced to an effective 
de-dimensional transport mo del in two steps (see Figure [1]). First, as we did in 
iBruna and ChapmanI (|2012bl) . the dimensions can be reduced from dN to d (indi- 
vidual to population-level description) by looking at the marginal density function 
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of one particle (the first particle, say) given by p{xi,t) = f P{x, t) dx2 • ■ ■ dx^r (the 
particle choice is unimportant since all the particles are identical). Second, we will 
exploit the geometry of the domain fi to further reduce the dimensionality by k, the 
number of confining dimensions. To this end, we will introduce the narrow- domain 
variables and obtain, from the d-dimensional density p(x, t) a reduced effective den- 
sity pe{xe,t), with Xe G R''^ • For cases (NC2) and (NC3), the effective density pe 
will be a one-dimensional density pe{x,t) along the channel axis. For (PP), it will 
be an effective two-dimensional density on the plane, pe = pe{x,y,t). 

For the sake of clarity we illustrate the derivation for the two-dimensional 
case (NC2) for both point and finite-size particles; the extension to the three- 
dimensional cases follows similarly and the respective models are only given in a 
summarised form. 



2.2 Point particles 

We begin by considering the case of point particles, for which the first reduction 
71 in Figure [T] from to one particle is straightforward. Since the particles are 
independent, P[x,t) = Ili^i PC-'^-i' ^'^'^ 

^(x,t) = Vx- [VxP-f(x)p] in (4a) 
0= [Vxp-f(x)p] -n on dQ, (4b) 

where n is the outward unit normal to dQ. Thus we move to the second model 
reduction T2 which is applied to Using the definition of fi pa|l . we want to 
exploit the smallness of H. We introduce a change of variables to the narrow- 
domain variables, which consist of rescaling by e the variables corresponding to the 
confined dimension: 

x = x, y = ey. (5) 

Introducing h such that H = eh, the domain i7 transforms into uj = [—5,5] x 
rescaled domain, we define p{i,t) = ep(x, f). (The factor of e is 
introduced so that both p and p integrate to one in their respective domains f2 
and ui.) Then ^ becomes 

''ft^^^'^ = ''§i (i - f'^^^^y^p) + 1(1 - 'y^p) ■ (6a) 

in UJ, with boundary conditions 

|| = /i(i,ey)p on S: = ±^, (6b) 

I? = e/2(i,ey)p on y = (6c) 
dy 2 

where /i and /2 are respectively the horizontal and vertical components of the 
external force f . Expanding p in powers of e, Taylor-expanding /i and /2 around 
{x,Q), and solving (|6ap with the boundary condition (|6cl) gives, at leading order. 
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that p is independent of y. Integrating (|6a|) over the channel's cross section and 
using (f6c)l we find that, to 0{e) 



dt 

r-h/ 



\x,t) 



d_ 

dx 



dpe 



- fi{x,0)pc 



fG [-1/2,1/2], 



(7) 



wliere pe = J_^y,pdy is the effective one-dimensional density along the channel. 
Tliis equation is complemented with no-flux boundary conditions at i; = ±1/2. 
Equation ^ can be generalised to three-dimensional geometries as 



dpe,. 

dt ^""^ 



Xe G OJe 



(8) 



with no-flux boundary conditions on cJoje, where Xe G a;e are the coordinates in 
the effective domain (i.e. the one-dimensional axis for (NC3) as in (I7|), or the two- 
dimensional plane for (PP)). The effective drift fe is the projection of the full drift 
vector onto the effective domain loe- The initial condition is pe(xe,0) = po(xe), 
where po(xe) = Jj-^iv -PQ(5)5(xe - xi,e)d£. 

A common extension to ([7} is to suppose that the channel has a non-c onstant 
cross section, h = h{x). The simplest model is the Pick-Jacobs equation ( Jacobsl 
Il967[) . which in our notation reads 



dpe 
dt 



{x,t) 



d_ 

dx 



h{x) 



d 



Pe 



dx \h{x) 



- fl{x,0)pe 



(9) 



and is valid for eh' {x) small. Generalisations to this equation to account for the 
chan nel curvature (a higher-order term) have been given in (jReguera and Rubj 
I2OOII) . The key step in deriving (|9]) is to assume that the full two- or three- 
dimensional probability density p(x, t) reaches equilibrium in the transverse di- 
rection, that is, it is assumed to factorise as p(x, f) ~ pe{x,t)p{'k), where p(x) is 
the local equilibrium distribution of y (and z, for d = 3), co nditional on a g iven x 
(the normalised Boltzmann-Gibbs probability density); see lZwanzid ( 1992h . 

In what follows, we keep h constant since the inclusion of a variable channel 
width in the analysis for finite-size particles is not straightforward. 



2.3 Finite-size particles 

We now describe the main result of this paper: the model of the effective dynamics 
in a confined domain for the dri ft-diffusion of finite-si z e parti cles. Using a similar 



technique to our previous work lBruna and ChapmanI (|2012bl) . we are able to re- 
duce the Fokker-Planck equation ([2]) for the joint probability density P{x, t) of A'' 
interacting finite-size particles in a confined domain H to the following effective 
equation for the marginal density pe{xe,t): 

^(Xe,t) = V^^ ■ { [1 + (iV - l)e'*=Q^Pe] V^^pe " /e(&e)pe} , (10) 

for Xe £ uJe C , where de are the effective dimensions of the reduced domain uje. 
The coefficient a/j, which depends on the geometry of the problem, determines how 
the excluded volume varies with the confinement parameter h. This equation is 
complemented with no-flux boundary conditions on duje and initial datape(a;e, 0) = 
Po{xe). Below we specify the coefficient for some specific cases. 
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1. Two-dimensional channel of width h (NC2), de = 1: 



/l2 



7r/i - - + ©(1 - /i) ( -(2 + - ft2 _ 2/iarccos(/i) 



(11) 



where ©(s) is the Heaviside step function, 



a; < 0, 
a; > 0. 



2. Three-dimensional channel of cross-section h x h (NC3), de = I: 



0(h - 1) ( -nh+j^]+Oil~ h)m{h) 



(12) 



where 



with 
s{h)-- 



i{h) = s{h) + < 



15 15 



0<h< 



71 <'^' 



(13) 



h[h'^ - Qh^ -h 4/i - 3) - 2/iarcsin(/i), 



-h -h^(h^ — 6)arccot f '^^^^ ^ _ — /i^arccot f ^ .^ ^ - 

6 ^ ' \ J 3 V2/iVl-2fe^ 



- arcco I i_5^2 + 6/i4 + 4/i2^(l-2/i2)(l-/i2) 



+ /iarcsin(/i) — ift(ft^ — 6/i^ — 3) arcsin 

3 — h'^ 



at{h) = — A/l-2fe2(/i^ + 9/1^ + 4) ^/i (^2/i^ - 12/1^ + 8/1-3 

^1,3. ,2 .^ ^ /^ 2feVl-2fe2 ^^ 4^2 , / 1 - 2fe^ - /i^ 
+ (/. -6)arccot^^-3^j+-ft arctan ^^^^^^=_ 

1, /4/iVl-2/i2(3;i2-l)\ 

+ 2'^^'^^'°* [ i-m2 + m4 j ■ 

3. Three-dimensional parallel plates a distance h apart (PP), de = 2: 

a/i = [(/i'(6 - /t')) ©(1 -h) + {8h - 3)0{h - 1)] . (14) 
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2 4 6 8 10 

h 

Figure 2 Excluded-volume coefiicient as a function of the confinement parameter h in 
three cases: two-dimensional channel (NC2) l lllll . three-dimensional channel (NC3) II12I I. and 
parallel plates (PP) llTill . 



The coefficient q/j corresponding to cases (NC2), (NC3), and (PP) is plotted in 
Figure [2l 

Equation (|10p describes the probability density for finding the first particle at 
position Xe at time t. Since the original system ([2]) is invariant to permutations of 
the particle labels, the marginal density function of any other particle is the same. 
Thus the probability distribution function for finding any particle at position Xe 
at time t is simply Npe- 

In Figure [3] we sketch the narrow-channel domain (rescaled by e) for various 
heights h. The physical domain (of width h -\- 1 \n the narrow-domain variables) 
is delimited by the solid black lines, while the configuration domain (of width h) 
corresponds to the yellow shaded region delimited by dot-dash lines. 

The nonlinear diffusion term in (|10p is proportional to the effective excluded 
volume created by the remaining (A*' — 1) particles as well as the domain walls after 
the dimensional reduction. For example, in the (NC2) case ea^ is the effective one- 
dimensional excluded interval, corresponding to the excluded area divided by the 
height of the cross section available to a particle centre (see Figure[3]). When /i = 0, 
a particle of diameter e excludes an interval of 2e (this explains why = 2 for 
ft = 0; see Figure[2]). As the channel width increases, the value tah decreases since 
the whole width of the channel is not always excluded by a given particle. As h 
gets large ea^ gives the ratio of the area excluded by the particle, yre^, to the cross 
section height eh, so that ~ Tr/ft as ft — >■ oo [see (|lip ]. 

While Q/j gives the effective excluded volume after dimensional reduction, the 
actual excluded volume is proportional to fta^ . This is plotted in Figure d] [along 
with the corresponding expressions for (NC3) and (PP)]. 

It is clear from Figure [3] that the excluded volume due to a particle varies 
depending on its position in the channel's cross-section: in (NC2), while a particle 
excludes an area of tt when it is far from the channel walls, it only excludes half 
of this area when in contact with the channel walls (less if ft < 1). This effect, 
known as an entropic effect, implies that the average excluded area over possible 
locations across the channel width decreases as the channel narrows (ft — >• 0). As 
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Figure 3 Sketch of the channel domain (shaded in yellow) for different values of h with par- 
ticles of diameter one (note that here we are depicting the actual particles, not their excluded- 
area which has radius one). Single- file channel for < h < 1 (with h = Q being the extreme 
case in which particles can only move in the axial direction). When h = I particles can just 
pass each other, and for h > 1 (bottom row) particles can more easily change order. 











(NC2) 


/ '' 


- - -(NC3) 




- - (PP) 



4 6 
h 



10 



Figure 4 Excluded volume ai^A as a function of the confinement parameter h in three cases: 
two-dimensional channel (NC2) (A = h), three-dimensional channel (NC3) (A = h^), and 
parallel plates (PP) (A = h). 



the channel width h grows, the boundary effects in which the excluded area is 
reduced contribute less and less to the average value, implying that the average 
excluded area tends to the constant value tt as /i — >• oo, which corresponds to the 
"bulk" excluded area. This is confirmed in Figure^ Similarly, tends to 47r/3 for 
the three dimensional cases as this is the rescaled excluded volume (the volume 
of the unit sphere). As /i — >• the average excluded area hafi — >• since in the 
extreme confinement cases almost all of the actual excluded area lies outside the 
domain available to a particle's centre. 
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2. 3. 1 Effective equation for the volume concentration 

In our derivation of (|10|) we do not require A'^ to be large: in fact, equation (|10|) 
is valid for any N (as long as the volume fraction is small), so that one could set 
= 1 or 2 if required. This equation gives the probability of finding a particle at a 
given position at a given time. However, for large A'' such that N — 1 ^ N we can 
introduce the volume concentration Ce = (j)pe, where is the total volume fraction 
of particles, and rewrite equation (|10p as an equation for the concentration of 
particles in the system^ 

^{Xe,t) = Vj,^ ■ [{\ + ghCe)V s,^Ce - f ^{xe) Ce] , (15) 

with 

(NC2): gh = -{h+lW, 

TT 

(NC3): gh = -{h+lfah, 
n 

(PP): g^, = -{h+l)aj„ 

TT 

These expressions for are plotted in Figure \5\ We note all three have a finite 
value as both — > and /i — >■ oo and, most importantly, that they have a relative 
maximum at h* > 1. In particular, h* = 1.47, h* = 1.28, and h* = 1.2 for (NC2), 
(NC3), and (PP) respectively. 

14 



12 



10 
6 



4 
2 

2 4 6 8 10 

h 

Figure 5 Coefficient a function of the confinement parameter h in three cases: two- 

dimensional channel (NC2), three-dimensional channel (NC3), and parallel plates (PP). 



The presence of this relative maximum at a fixed volume fraction is interesting, 
as it implies an optimal ratio between the particles' size and the confinement 



iVTre 
4(ft + l)' 

Nne 
6(ft-hl)2' 

iVTre^ 
6(ft + l)' 



(16a) 
(16b) 
(16c) 





(NC2) 




- - -(NC3) 




- - (PP) - 


/ ^ 

/ ^ 

- 'l ...... 

fi- 
ll : : 

[l 

1 

1 







^ Note the factors (1 -|- h) in <j>: this is because </) is the total volume of particles divided by 
the actual volume of the channel, not the volume available to a particle's centre. 
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dimension at which excluded-volume effects, and thus the effective transport, are 
maximised. In terms of the physical domains, it corresponds to a narrow domain 
around 2.2 to 2.5 times wider than the particles' diameter, so that two particles 
can just diffuse side by side in the channel. 




Figure 6 Sketch of channel of width 5e with and without an intermediate wall that creates 
two lanes. The domain available to the particles centres' in each case is shaded in yellow. Left: 
one narrow channel with h = 4. Right: two narrow channels each with h = 1.5 (roughly equal 
to h* for maximal exclusion effects). 



Thus the theory predicts that diffusive transport in a (NC2) channel of width 
5e {h = 4) may be increased by dividing the channel into two sub-channels of 
width 2.5e {h = 1.5), as shown in Figure[6l This operation gives an increase in gf^ 
of 7%. The increase is more dramatic in the three-dimensional case: if a (NC3) 
square channel has an original width and depth of 4.6e (so that h = 3.6), then 
subdividing it into four identical channels of width 2.3e (so that each has h = 1.3) 
gives a relative increase in is of 19%. 



3 Limiting cases: from single-file diffusion to unconfined diffusion 

We have briefly discussed the limiting behaviour of as h and h ^ oo above. 
Here we examine these limits in equation (|10|1 . and check that they agree with 
existing results. For a channel of width h (NC2 or NC3) equation (|10p interpolates 
between two limiting cases: a single-file channel {h — > 0) and an unconfined two- 
(NC2) or three-dimensional (NC3) domain {h — > cxd). For the (PP) case, the limit 
h gives an (unconfined) two-dimensional diffusion, while the limit h ^ oo 
gives three-dimensional unconfined diffusion. Finally, the extension of the square 
cross section (NC3) to a rectangular cross section hxm can be used to interpolate 
between (NC2) (as m — >■ 0) and (PP) (as m — > cxj). 

In this section we will examine these limits by comparing the limiting behaviour 
of our model (|10p with the limiting problem of diffusion of hard spheres in R** 
for d = 1,2,3. This problem has been studied extensiv ely, especially in the one- 
dime nsional case, which is known as single- file diffusion ( Lizana and AmbiornssorJ 



l2009f ). For the cases d = 2,3 we wil l use the results from our p revious work of 
unconfined diffusion of hard spheres ( Bruna and Chapman|[2012lJ ). 
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3.1 Limit to an unconfined domain: h ~> cxd 

As h increases we expect the bound ary effects contained i n to vanish and to 



recover the "bulk" equation found in lBruna and ChapniianI (|2012b[) : 



^(x,t) = Vx- {[f + (iV-l)ae'^p]VxP-f(x)p} in O, (17) 

where i7 C M'* as given in Q, a = n for d = 2 and a = Att/S for d = 3 (or 
a = 2{d — l)Tv/d for d = 2,3). It is important to note that this equation is only 
valid for H = 0(1) (that is, when !? has volume order one). 

In order to take the limit /i — >• oo in our model (|10p . it is convenient to use the 
original density p in rather than the effective density pe. In other words, we 
consider the following equation for p = pe/A (where A is the cross-sectional area): 

^{xe,t) = Va^-[[i + {N- l)a!,Ae''^p\ VAj-f,{xe)p] , (18) 

where af^A is the excluded- volume coefficient shown in Figure [4] As can be seen 
in the figure [or in the formulas for (|lip - p4p ]. the limit of Of^A as — cx) is vr 
for the two-dimensional channel (NC2) and Att/S for the three-dimensional cases 
(NC3) and (PP). Therefore, the limiting behaviour of psp as /i — > oo corresponds 
to replacing af^A by a, where the latter is given in the bulk equation p7|) . The 
last step to show that p7|) is indeed the limiting model of psp is to integrate (|17p 
over the cross section to reduce it to a de— dimensional equation as (|18p . In other 
words, the limit /i — ^ oo of the confined-domain model psp should coincide with 
the limit _ff — >• of the bulk equation p7p . Rescaling the confined dimensions by e 
[cf. ([5])] and integrating (|17p over the cross-section, it is straightforward to arrive 
at the following equation for p = e^p: 

^(ie, t) = V^^ ■ { [l + (iV - l)ae'''p\ V^J - /,(ie) p} , (19) 
where we have used no-flux boundary conditions on the cross-section boundaries. 



3.2 Limit to single-file diffusion: — > 



We now consider the limiting case /i — >• in the narrow channel cases (NC2) and 
(NC3). From the plot of in Figure[2]we have that lim/j_j,g a/j = 2 in both cases. 
This is confirmed by taking the limit /i ^ in ^ for (NC2) or ^ for (NC3). 
Thus the single-file limit of (|10p is 



dpe 



[l + 2{N~l)epe]^~hix)p 



(20) 



We see that the effective diffusion coefficient for large (such that — 1 ^ 
A'') is -D'^(c) = (1 -|- 2c) with c = Ntpe being the particle concen tration, which 
is consistent with that derived by lAckerson and Fleishmai] (|l982l) for a uniform 
particle concentration c = Ne. 

We now compare (|20p with the classic one-dimensional model of diffusing hard 
rods. It is well known that the one-dimensional diffusion of finite-size particles 



14 



Maria Bruna, S. Jonathan Chapman 



can be mapped onto a point -particle problem (cf . iLizana and Ambi6rnssonll2009l) . 
Using this trick, a fast diffusion equation for the evolution for the marginal density 
of A'' rods of length e under no external force (/i = 0) a nd in t he the rmodynamic 
lim it (A^ — >■ (X), L — >■ oo, N/L 4) finite) is found in iRostI (|l984l) ( in French, 



see iBodnar and Velazaue j l2005l for an explanation of Rost method in English) , 
namely 



dt ' Si Y (1 — ep)-^ dx 

where p = Npe is the number density. Expanding the equation above in e we 
obtain, to 0{e), 

f(^,.) = ^((l + 2iV.^.)f), (22) 

which is in agreement with the large A'' limit of (|20p . An alternative derivation of 
(|20p using matched asymptotics on t he original pro blem (without eliminati on of 
the ha rd-core parts) can be found in iBrunal (|2012l) . It differs from that in iRostI 



( 1984) in that it is valid for any N and allows an external force field f. 



3.3 Other limits 



From (|14p we see that lim^_^Q = tt (see also Figure [2|), from which it is straight- 
forward to show that the limit ft — >■ of (PP) corresponds to an unconfined two- 
dimensional diffusion. 

A final limit to consider concerns the generalisation of the three-dimensional 
channel (NC3) to a rectangular cross section h x m. We have already seen above 
that this tends to the single-file diffusion model for ft = m — >• 0, and to an uncon- 
fined three-dimensional diffusion for ft,m — >• oo. Now, keeping ft fixed, the extra 
parameter m will allow us to interpolate between a two-dimensional narrow chan- 
nel (NC2) of width ft as m — s> and two parallel plates a distance ft apart (PP) 
as m — ^ oo. For m > 1 it can be shown that the rectangular counterpart of (|12|) 
reads: 



1 
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+ e(ft- 1) 



4n 



hm - 



-(ft -I- m) 



+ 0{l~h)s{h,m) \ , (23) 



with 



s(ft, m) = Trmft j 1 



ft arcsin ft + 
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(2ft'' - 9ft^ 



Now, in other to compare between the one-dimensional (NC3) model and the 
two-dimensional (PP) model, the relevant quantity is ma/j^ (so that the one- 
dimensional effective density pe is mapped onto a two-dimensional plate of width 
m). We find that 



,. (NC3) (PP) 

hm ma, = a, ' 
hm h ' 



as expected. 
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4 Model analysis 



4.1 Numerical analysis of time dependent solutions 



In this section we compare solutions of our effective nonlinear diffusion equation for 
the (NC2) case with direct stochastic simulations, averaging over 2 • 10^ stochastic 
realisations of the corresponding individual-based model ([1}. For the PDE, we use 
the method of lines with a standard second-order finite-difference discretisation of 
the spatial derivatives. For the coupled system of SDEs, we perform Monte Carlo 
(MC) simulations of the discretised version of ^ using the Euler-Maruyama 
method, 

X,{t + At)=X,{t) + f{X,{t))At + V2Ati„ (24) 
where ^^ is a 2-vector whose entries are independent normally distributed random 
variables with zero me an and unit varian ce. The reflective boundary conditions on 
dn implemented as in lErban et all ( 20071) . namely, the distance that a particle has 
travelled outside the domain is reflected back into the domain. Care must be taken 
for very narrow channels to account for the possibility that a particle has travelled 
outside the domain by more than a width h. The particle-particle overlaps are 
implemented similarly: the difference e— ||Xi(t-|- At) — Xj(t + At)\\ corresponds to 
the distance that particles have penetrated each other illegally. Then we suppose 
that each particle has travelled the same illegal distance, and we separate them 
accordingly along the line joining the two particles' centres. This approach works 
well for low volume fractions, but may run into difhculty at high volume fractions 
when the separated particles may suffer further overlaps. In that case the algorithm 
iu iScala et aL (|2007i) . an event-driven Brownian dynamics, becomes more suitable. 




Figure 7 One-dimensional density pe{x,t) in a channel of width eh with no-flux boundary 
conditions at a; = ±0.5; solution of the continuum model l llOII for finite-size particles (solid blue 
line) versus individual-based model simulations (blue circles) and limiting continuum models 
(point particles jTj, cross-solid green lino; single-file limit II20I I. dash red line; unconfined limit 
I I19II . dot-dash black line), (a) Initial (f = 0, cyan) and final (t = 0.05, blue) densities pe (lines) 
and histograms (data points), (b) Final density and histogram, together with limiting cases. 
Parameters are h = 3, e = 0.01, and TV = 30, At = 10'^ . 



Figures [3 and [8] show the numerical results at t = 0.05 for h = 3, e = 0.01, 
A'^ = 30, f = with no-flux and periodic boundary conditions at the channel ends. 



16 



Maria Bruna, S. Jonathan Chapman 



respectively. At initial time, the particles are uniformly distributed in a segment of 
length 0.2 [Figures[3a) and[8l[a) in cyan triangles and dash line]. The data points 
show the one-dimensional histogram obtained by averaging the MC results over the 
channel's cross section. To test the importance of the excluded-volume interactions 
and the confinement, we also compare with the corresponding solutions with point 
particles and the limiting cases as — ;> oo and — > of Subsections 13.11 and 13.21 
respectively. 



(a) 


1 
1 
i 

1 








Figure 8 One-dimensional density pe{x,t) in a channel of width eh with periodic boundary 
conditions at a; = ±0.5: solution of the continuum model lllOII for finite-size particles (solid blue 
line) versus individual-based model simulations (blue circles) and limiting continuum models 
(point particles lO, cross-solid green line; single-file limit Il20ll .dash red line; unconfined limit 
I I19I1 . dot-dash black line), (a) Initial {t = 0, cyan) and final (t = 0.05, blue) densities pe (lines) 
and histograms (data points), (b) Final density and histogram, together with limiting cases. 
Parameters are /i = 3, e = 0.01, and N = 30, At = 10"^. 



In both cases we see very good agreement between the stochastic average and 
the solution of the narrow-channel model pe, whilst there are noticeable differences 
between the three limiting models, namely the point particles, single-file and bulk 
limits. In order to quantify the error committed by the limiting models, we note 
that, for the chosen parameters, the volume fraction is ^ 6% and the nonlin- 
ear coefficient in p6ap is g/j 4.6. The corresponding nonlinear coefficient in the 
limiting cases is for the point particles limit, 2 for the single-file limit and 4 
for the bulk limit. The different between the narrow channel coefficient and these 
limiting values are consistent with the differences observed in the numerical solu- 
tions: while for /i = 3 the bulk limit compares reasonably well with the stochastic 
simulations, the single-file and the point particles limits show more important dif- 
ferences. However, we note that our model gives the best agreement with the MC 
simulations results. 

The results of Figures [7] and [8] suggest that, depending on the value of h, either 
the single-file or the bulk limits will be more appropriate, while our narrow-channel 
model captures the whole range of h from to oo. To investigate this further, in 
Figure [9] we examine how the narrow-channel model compares with three limiting 
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cases (single-file, bulk, and point particles) for various values of the channel width 
h while keeping the volume fraction 4> fixed0 



ft = 0.5 h = l 




X X 

Figure 9 One-dimensional density pe{x,t) in a two-dimensional channel of width e/i with 
fixed volume fraction (p = 0.05 at time t = 0.05. No-fiux boundary conditions at x = ±0.5 
and uniformly initial conditions in \x\ < 0.1. Solution pe of the narrow-channel equation lllOI I 
(solid blue), versus the point particles limit jTj (cross-solid green line), the single-file limit II20I I 
(dash red line), and the bulk limit II19II (dot-dash black line, shown curve is hp). 



As expected, for h = 0.5 the narrow-channel solution agrees very well with the 
single- file limit, since in this case the cross-sectional space is not enough to let 
particles pass each other (see top right plot in Figure[3l). In contrast, the bulk case 
solution is far apart from the previous two, since the approximation that boundary 
effects are negligible is poor for h = 0.5. As we increase h, we can observe how the 
single-file solution (in dash red) moves apart from the narrow-channel solution (in 
solid blue) , while the bulk case solution (in dot-dash black) becomes closer to the 
latter. When h = 5 the narrow channel and bulk curves are nearly overlapping 



^ We keep the volume fraction <f> = Nne/4:{h + 1) fixed by varying e as h changes. 
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each other, indicating that at this channel width the boundary effects are almost 
negligible. These observations can also be made by looking at the graph of in 
Figure [SI in particular by considering the difference between at a given value 
of h and the extreme values at ft = 0, go = 2, and at h 



4.2 Stationary solutions in periodic channels 

In this section we study the steady states of our continuum model and use the 
narrow-channel equation as an extension to the one-dimensional linear ratchet 
model for Brownian motors. For convenience, we focus on the one-dimensional 
equation (|15p in terms of the coefficient g/^ for the narrow channel cases (NC2) 
and (NC3) only. As in the time-dependent case, the numerical solutions of the 
PDE are compared with corresponding stochastic simulations of the individual- 
based model ([1]). We suppose that the force along the channel axis fi{x) is the 
gradient of a potential V{x), so that fi{x) = —V'{x), where the prime indicates 
differentiation. Then (1151) can be written as 



~dr ~^ di^^''^^ ^ where u = - — {\ogpe + ghipPe + V{S:)) . (25) 
Th e quantity u can be interpreted as a flow dow n the gradient of the free energy 



J" ( Carrillo et all [20031) associated with psp : see iBruiial ( 20121) for more details 



The stationary solution of (|25p with no-flux boundary conditions is obtained by 
minimising the free energy, which corresponds to solving 

logpe + ghfpPe + V{x) = C, 

with the constant C determined by the normalisation condition on pe- For the 
application to ratchet systems, we are interested in periodic solutions of (|25p with 
a (constant) flux Jo = peiilfl 

The one-dimensional Fokker-Planck equation dT]) for point particles with /i {x) = 
— V'{x) can be used to model directed particle transport under ratch et potentials 
(i.e. potentials spatially asymmetric with respect to their maxima I Slater et all 
Il997l) . These potentials may describe a periodic asymmet ric free-energy substrate 
in th e case for molecular motors through microtubules ( Kolomeiskv and Fished 



l2007l) or steric interactions in the case of polyelectrolytes (|Slater et allll997l ) In 
particular, theoretical approaches to molecular motors such as kinesin have fo- 
cused on either a one-dimensional continuum model such as ([7]) (thus ignoring the 
interactions between different motors and the other dimensions) or on stochastic 
simple-exclusion models on a lattice ([Kolomeiskv and Fished l2007l) . More gener- 



ally, in these applications the modelling has focused on the form of the ratchet 
potential in order to make the model more realistic, instead on the possible inter- 
actions between the particles involved in the transport. In the remainder of this 
section we will explore the effects that interactions (speciflcally, excluded-volume 
interactions, but this could be extended to other types of interactions) can have 
on such models. To this end, we consider a specific ratchet model with a ratchet 



For periodic boundary data, Jo is an extra degree of freedom determined by imposing 
periodicity. 
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potential consisting of a periodic part plus a constant external force or Ult. We use 
the tilted Smoluchowski-Feynman potential 

V{£, Fo) = sin(27r£) + 0.25 sin(47r£) - Fqx. (26) 

Two plots of V{x) for different values of the tilt Fq are shown in Figure [TOl It 
can be shown that in the long-time limit the sign of the particle current (or net 
motion) agrees with the sign of the tilt F g and that, in the absence of tilt (Fq = 0) 
there is no net motion ( Reimannl |2002|) . An interesting feature of this model is 



that the relationship between the tilt Fq and the flux Jq is nonlinear (see thick 
black line in Figure [T2|l . 




Figure 10 Tilted Smoluchowski-Feynman ratchet potential V(x,Fo) in Eq. II26II for Fq = —I 
(left) and Fq = -3 (right). 



Next we examine how the Fq vs. Jq relationship changes when nonlinear effects 
(due to the finite-size of particles and confinement) are included in the model. To 
this end, we solve for the stationary states of (|25p when V{x) is given by (|26p for 
various tilts Fq and compute the stationary flux Jq of the resulting solution. For 
each tilt Fq, we must find the flux Jo and the stationary solution pe{S:) such that 
the periodic and normalisation conditions are fulflUed: 

{^+ 9h<l>Pe)pe + V{x,Fo)pe = -Jo, 

Pe(-l/2)=Pe(l/2), (27) 

J-i%Pe dx = l. 

We solve this problem numerically using Chebfun ( Trefethen and others! [20 



in MATLAB. Solutions of (|27p for an increasing value of the constant force Fq 
(varying from —6 to +6) are shown in Figure 1111 The left panel corresponds to 
point particles (g/j = 0), while the right panel corresponds to finite-size particles 
with giicj) = 0.6. The diagram of the resulting steady fiux Jq versus the tilt Fq is 
shown in Figure [12] for increasing values of gfi4>. We observe that the relationship 
is nonlinear for point particles (g/j = 0, thick black line), but it appears to become 
linear as excluded- volume effects get larger (i.e. as g/j^ increases). This is physically 
reasonable, since point particles get easily trapped in the wells of the potential, 
even if these wells are relatively small. In contrast, it is easier for finite-size particles 
to escape, as they may not all fit in the potential well and the nonlinear diffusion 
makes the potential barrier easier to overcome. 
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Figure 11 Stationary solution Pe{x) under the tilted Smoluchowski-Feynman potential II26I I 
for various values of Fq. Stationary solutions of II25I I for point particles, ffh = (left) and for 
finite-size particles with g;,</i = 0.6 (right). 




Figure 12 Steady-state flux Jo versus tilt Fq from solving II25II for increasing values of gfi<p 
(from to 1). The coloured lines correspond to finite-si ze particles (g^ > 0). The thick black 
line corresponds to point particles {gf^ = 0) as shown in iReiman 3 1I2OO2I . Figure 2.4). 



In Figure [T51 we compare the stationary solutions pe for point particles (g^ = 0) 
and finite-size particles (with an excluded-volume coefficient of gi^ = 0.15) in two 
tilting scenarios: for Fq = —6 and for Fq = 2.5. We observe that while the solutions 
are almost overlapping for a tilt of _fo = —6, they are considerably different for 
Fq = 2.5. This is because for Fq = —6 the potential V is so tilted that it ceases to 
have a local minimum within each period. As a result, the "advantage" of finite- 
size particles that could more easily overcome the local minima in the potential is 
lost. 

We conclude this section by comparing the equilibrium solutions of the contin- 
uum model with the corresponding stochas tic simulations of the discr ete model. 
We use the Metropolis-Hastings algorithm ( Chib and Greenberg|[l995l ) to sample 
from the stationary density of the full-particle system, and compare the resulting 
histogram (averaged over the cross section) with the stationary solutions pe of 
(|25p . We consider the (NC2) case for which the coefficient gii is maximised for a 
fixed volume fraction that is, a channel of width h* = 1.47. The other param- 
eters used in the simulations are e = 10"'^, A'^ = 133 and Fq = 2.5. Figure [Til 
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Figure 14 Histogram of the stationary density p(x) under the potential V{x, Fq) in II26II for 
point particles (top plot) and finite-size particles (bottom plot). Parameters are Fq = 2.5, 
h = 1.47, e = 10"^, N = 133. Histograms computed by 10'' MH steps. 



displays the histograms after 10^ steps of the MH algorithm. The histogram for 
point particles (upper plot in Figure [Ti|l does not vary in the cross-sectional direc- 
tion, as expected. In contrast, the histogram for finite-size particles does display 
some variation in the y-direction: more particles want to be near the boundaries 
than in the centre of the channel. This is because a hard-disk particle on the 
boundary excludes only half of the area that would exclude in the centre of the 
domain (recall that the channel of width h is the domain available to the particles 
centres, not the physical domain), and this is entropically favourable. 

We plot the theoretical predictions by solving (|25p for both point and finite- 
size particles alongside their respective simulation counterparts, and observe an 
excellent agreement in both cases. We examine the importance of taking into 
account the actual width of the channel, which in this example is only h = 1.47, 
by solving the analogous stationary problem for the single-file model (|20p . We plot 
the result (dot-dash black line in Figure [T5)) and observe that the single-file model 
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overestimates the excluded-volume effects, as demonstrated by the ffatter density 
profiie. 



3.5 - * ^ 

/ \ 
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Figure 15 One-dimensional stationary density p<i{x) under the potential _Fo) in I I26II . 

Solutions of II25II for point particles (g^j = 0, dash red line) and finite-size particles (g^ = 0.1, 
solid blue line), with Fq = 2.5, N = 133, h = 1.47, e = 10"'^. Cross-sectional averages of the 
histograms in Figure [T4l for point particles (red asterisks) and finite-size particles (blue circles). 
Stationary solution of the single-file equation II20II is shown in a dot-dash black line. 



5 Derivation of the narrow-channel equation (llOp for the (NC2) case 

This section is devoted to the derivation of (|10p in the two-dimensional channel 
(NC2) case. The derivation of the three-dimensional cases (NC3) and (PP) or 
other (simple) geometries follows similarly (see Subsection 15.31 for an outline of 
the conditions/calculations to be carried out). 



5.1 Step A: Reduction of the number of particles, from individual- to 
population-level 

Our starting point is the Fokker-Planck equation for TV hard-disc particles (|2ap . 
defined in the high-dimensional (configuration) space C K^^. Recall that the 
configuration space has holes that correspond to illegal configurations (with par- 
ticles' overlaps) on which the no-flux boundary conditions (|2b[) hold. The aim 
of this subsection is to derive the corresponding Fokker-Planck equation for the 
one-particle density p(xi,f) = J P(i, t)dx2 . . .x^v, where P{x,t) is the joint prob- 
ability density of the A*' particles. We first note that the integration of ^ over 
X2 , . . . , x jv results in int egrals over contact surfaces on which P must be evalu- 



ated ( Bruna and Chapman 2012b}). When the particle volume is small, the dom- 



inant contributions to these contact integrals correspond to two-particle interac- 
tions, so that we can set N = 2 and then extend the result to A'' arbitrary in a 
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straightforward manner. How ever, in contrast with the unconfined case studied in 
iBruna and ChapmanI (j2012bl) , under confinement conditions the particle-particle- 
wall interactions (three-body) are not negligible and must be taken into account. 
For two particles at positions xi and X2, ([2]) reads 

-^(xi,X2,t) = Vxi-[VxiP-f(xi)P]-hVx,-[Vx,P-f(x2)P] in «e', (28a) 

0= [VxiP-f(xi)P]-ni + [Vx,P-f(x2)P]-n2, (28b) 

on Xi G dil and ||xi — X2II = e. Here = n,;/||ni||, where nj is the component of 
the normal vector n corresponding to the ith particle, n = (ni, 112) ■ We note that 
ni = on X2 G dO, and that fii = —112 on ||xi — X2II = e. 

5.1.1 Step A.l: From N particles to 1 particle 

We denote by f2(xi) = Q \ Be(xi) the region available to the centre of parti- 
cle 2 when particle 1 is at xi. Note that when the distance between xi and 
dQ is less than e the area |i7(xi)| increases (see Figure I16|) because the area 
Z//(xi) = Pe(xi) n Q excluded by particle 1 changes with xi. The points on which 
the two particles are in contact are given by the collision boundary 

Cxi = {X2 G 12 s.t. ||X2 -xill = e} . (29) 

Integrating Eq. (|28a|) over i7(xi) using the Reynolds transport theorem (on 
the moving boundary Cxi), the divergence theorem, and the boundary condition 
([28bll yields 



^(xi,t) = Vxi • [VxiP-f(xi)p]-F^ -(VxiP + Vx,P)-n2d52. (30) 



We denote the collision integral above by T. If we now consider the case of A*' 
particles we would obtain a collision integral for each pair, so that after some 
particle relabelling the corresponding equation is 

^(xi,t) = Vxi-[VxiP-f(xi)p] + (iV-l)/" -(Vx,P + Vx.P)-n2dS2. (31) 

>^ Cxi 



eh 



_Ie/2 

Xl^ 



- 2e 



Figure 16 Sketch of the original channel geometry (solid black lines) and the eff'ective con- 
figuration space for the centre of a second circular particle, given by the boundary function 
9J7(xi) (dash red lines) which depends on the position of the first particle xi. The collision 
boundary function Cx^ forms part of the effective boundary 9J7(xi). 
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Equation (|3ip is halfway through transformation 71 (cf. Figure [T|), since the 
first half of the equation depends only on xi while the integral X still depends on 
the two-particle density P near the collision surface Cxi ■ 

At this stage, it is common to use a closure approximation su ch as P(xi, X2, t) = 
p(xi, f)p(x2, to evaluate! and obtained a closed equation forp ( Rubinstein and Kelle j 
Il989l) . However, the pairwise particle interaction — and therefore the correlation be- 
tween their positions — is exactly localised near the collision surface Cxi ■ Instead, in 
the next section we will use an alternati ve method based on matched asymptotic 
expansions to evaluate X systematically (jBruna and Chapmanir2012bh . 



5.1.2 Step A. 2: Matched asymptotic expansions 

We suppose that when two particles are far apart {\xi — 12] 3> 1) they are indepen- 
dent (at leading order), whereas when they are close to each other {\xi — X2\ ~ e) 
they are correlated. We denote these two regions of the configuration space fi^ 
the outer region and the inner region, respectively. We use the s-coordinate to 
distinguish between the two regions because the inner region spans the channel's 
cross section. Importantly this implies that the outer region is disconnected. 



Outer region In the outer region we consider the change to the narrow-domain 
variables ([5]) and define P(xi,X2,i) = e^J'Cxi, X2, t). This scaling is consistent 
with that introduced for the one-particle density in Section [2. 2 1 and is such that P 
and P each integrate to one in their respective domains and the narrow-domain 
variable equivalent la^. Then (|28a|l becomes 



(ii,X2,t) = 4(||-e/2(il,.yi)P)+4(£-e/2(i2,ej/2)P) 

for (xi,X2) G i^e, with 

h{^^,ey,)P = G 
§P--ef2{x„m)P = 



(32b) 



on 



-I- r, , 



(32c) 



for i = 1,2. The boundary condition on the collision line Cxi disappears for 
j£i — X2\ > e, so that it is "invisible" to the outer region|3 

In the outer region, we define Pout(xi, X2, f) = P(xi,X2,f) and look for an 
asymptotic solution to (I32|) by expanding Pout in powers of e. We find at leading 
order that Pout must be independent of the vertical coordinates yi and y2. By 
independence in the outer region we suppose that the leading-order solution is 
separable, so that it is of the form q{xi,t)q{x2,t) for some function q. Solving for 



* To see this, we write I I29I I in terms of the narrow variables. It becomes Cxj^ = {x2 6 ui 
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the next two orders in e we find tliat tlie solution in tire outer region is, to 0{€^), 

Pouti^i, X2,t) = q{xi,t)q{x2,t) 

+ e {17(^1, t)qix2, t) [2/1/2(^1, 0) + ^2/2(^2, 0)] + Ti[xi,X2, t)) 

+ ^'{5 [df^(^i,0)+y|^(^2,0)+(yi/2(ii,0)+y2/2(i2,0))'] (^3) 

+ Tl{xi,X2,t) [^1/2(^1,0) + 2/2/2(^2,0)] + T2(il,X2,t)|, 

where the Ti{x\,X2,t) are arbitrary functions of integration, determined by solv- 
ability conditions at higher order. Note that the invariance of P with respect to 
a switch of particle labels means that in the outer region both particles have the 
same density function q. From the solvability condition on the O(e^) terms above 
we obtain the following equation for q: 

Inner region In the inner region we introduce the inner variables 

xi = xi, yi = eyi, X2 = Xi + ex, 2/2 = ey2, (35) 

and define P(xi,X2,i) = e^f'(xi, X2, t). The contact boundary Cxi (|29|1 becomes 

Cy, = {{i,y2)eRx[~h/2,h/2] s.t. x'' + {y2~yif = 1], (36) 

and problem (|28p is transformed to 

e^f (xr,i2,t) = + + ^ - - [f2ix, + ex,ey2)P] 

+ [/i (^ 1. ^yi) - /i + fya)] P} (37a) 

with 

2i|| + (y2-2/i)(||-£) =.i||+.i[/i(ii+.£,ey2)-/i(ii,.2/i)]P ^^^^^ 

+ e(y2 - J/i)[/2(il+ei, £^2) - f2{xi,eyi)]P, 

on and 

9P- = ef2{ii,eyi)P on = ±|, (37c) 

§t = ef2{ii + ex,ey2)P on ^2 = ±|. (37d) 

In addition to (j37b|) - (|37d|) . the inner solution must match with the outer as 
X — >■ ±00. Expanding the outer solution in terms of the inner variables, which 
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corresponds to replacing xi = xi, yi = yi, £2 = xi + ex and 1/2 = y2 in (f33ll . and 
subsequently expanding in powers of e, we obtain the following matching condition: 

P^q^+e [xq-^ + (yi + y2)f2q^{il) + ri(ii,ii)] 

+ 6^ [f 90 + i{yi + y2)f2qqi, + iV2^q^ + 

/ , ' (37e) 

+ \ [{yl+yl)i^ + {Vi+V2ffl) q^ + {yi+V2)hTi{ii,ii) 

+ 5^2(2:1, +•• • as i— >-±cxj, 

where q = q{xi,t). Expanding P in powers of e, P ~ + eP^^^ + e^pCa) ^ ^ 

we find that the leading- and first-order solutions of (|37p are 



P^°^=q\il), (38) 
p(i) = i5(ii)^(ii) + (yi -Fy2)/2(ii,0)g^(5i) ri(ii,ii). (39) 

Unlike in the bulk problem (jBruna and Chapman|[20f2allbl ). the narrow-channel 
problem requires computing the second-order inner solution. The solution proce- 
dure is rather cumbersome and is omitted here. It involves a further change of 
variable x = y/2s to turn the problem into a Poisson problem, and solving two 
sub-problems numerically using the commercial finite-element solver CO MSOL 
Multiphysics 4.3. The second-order solution of (|37p is (see Appendix C.l in lBrunal 
I2OI2I for fuU details) 



'(2) = + l{jjl+yl)q^(^^+ /I) + y,y,f^q^ 

+ ^ {y^hiSk + y^^'-^) + {^-^) <i'Q2{i,yi,y2) (40) 



+ V2\q^ - {-^y - ^q^\Q^{x,n,y2) +T2{£i,x^) 



with Qi{x,yi,y2) = Vi{x/V2,yi,y2), where ui (s, yi , 1/2) and 1)2(5,1/1,^2) are given 
by (numerical solutions of) 



(41) 



~ 2 

V «i 


= 0, 






Viii • 


-2 
= S 


on 






= 


on 


= ±2 ' 


Vl 


-DilSl 


as 


s — >■ ±00, 



and 



~ 2. 

V V2 =0, 

VS2 ■ i> = 5(yi -m), on Vy^ , (42) 

Vi32-J> = 0, on = ±^, 

C2 ~ 0, as s — >■ ±00. 
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Here V stands for the gradient operator with respect to the position vector 
(s, yi,y2), 'Dy-^ is the transformed coUision boundary Cy-^ (I36p . and v is the outward 
unit normal on this mapped boundary, v = — -^(2s, yi — j/2, j/2 — j/i)- Finally, the 
constant field Di at infinity for vi is related to the integration function from the 
outer Ti, 

lim^.^^j ||i(ii,f2) 

D\ — 



a^q / dq \ 2 dfi 2 



Out of the derivation we also find that Ti{xi,x2) = Ti(|2;i — X2\) with Ti(x) 
differentiable satisfying Ti(0) = 0. In contrast, the contribution from the other 
outer function T2 is left unknown (it could be determined by matching higher 
order terms), but we are able to ignore it as it has a zero contribution to the 
collision integral I (as we will see in the next section). We note that, in (|40|) . q, 
fi and /2 are functions of the "outer" variable xi only, namely, q = q{xi,t) and 
= /.(ii,0). 

Combining (fSS]) . (|39)) and (f30|) we have the solution to the inner problem (f37)l 
up to O(e^). 



5.1.3 Step A. 3: Collision integral 

Now we go back to Eq. (|3ip and use the asymptotic solution of the previous sub- 
section to turn it into an equation for p(xi, t) only thus completing transformation 
71. Note that, since the integral I is over the collision boundary Cxi , it lives in the 
inner region and we must use P to evaluate it. 
In terms of the inner variables, T is 




where Cy-^ is given in (|36p and dl is the line integral along this curve (for yi fixed, 
see Figure [T7)l . Depending on the channel width h (relative to one, which is the 
radius of Cg^), the integration is over the whole circle or a part of it. Introducing 
the distances l\ = max(— /i/2 — yi,—V) and I2 = min(/i/2 — t/i,l), the angles at 
contact with the lower and upper channel walls are 6\ = arcsinZi and 62 = arcsin/2, 
respectively. (These are equal to ±7r/2 for no contact.) 

Writing 1 = e-2(x(°) + el^^'> + e^l^'^'> + • • • ) and using (gS]), we find that 

2:(°) = 0, (44) 

1^^^ =2f2q^^I0{h,yl), (45) 
1(2) = 2q^ (^+2/1) yiMh,yi) + (2g|^ - g2^) f^i{h,y^) 

+ + 2/1) M2(fe,yi) + 2 - ^) q'j[Q2]{h,y,) (46) 

+ '^^[l^,~{Si-f-^l']jm(.h,y,), 



28 



Maria Bruna, S. Jonathan Chapman 





- ~ ^ ^ 












k j 





h 



Figure 17 Domain of integration Cy^ (solid green line), distances l\ and Z2 between the 
vertical coordinate of the first particle and the lower and upper channel walls and corresponding 
angles Q\ and 92- 



where = /i(a;i,0), 

/xo(/i,yi)= / (y2-yi)d[=2(vT^- VT~i), (47a) 

l^i{h,yi) = / dl = l2\/l — I2 — ^ if + arcsini2 — arcsin/i, (47b) 

IJ-2(h,yi) = / (2/2 - yi)^ dl = ^ia/1 - If - h\/l ~ arcsin/i + arcsin/2, (47c) 
JCy, 

and J is the integral operator J[Q]{h, yi) = [<3g2(j'2 — yi) + Qxx] dL The terms 

J7[Qi] and J '\Q2\ are evaluated numerically with COMSOL (see Appendix C.2 in 
lBrunall2012l for more details). 

5.1.4 Population-level Fokker-Planck equation 

Combining (|45p and (|46|1 we obtain the first two terms of the asymptotic expansion 
for I, which depends on botlr tire channel width h and the elevation of the first 
particle yi but is independent of the position of the second particle. Thus we can 
drop the first particle label (the subindex 1) for clarity of notation. Inserting this 
expansion into (|31|) . we obtain an equation for the first particle 

^(x,f) = Vx-[VxP-f(x)p] + (7V-l)(e-il(i)+x(=^)) in n, (48) 

which involves the marginal density p(x, t), the outer density q{x, t) and the channel 
width h. This concludes the transformation 71 from N particles to one particle (see 
Figure [1]). 

5.2 Step B: Reduction of the number of geometric dimensions 

Following a similar procedure to that for point particles in Section 12.21 we will 
reduce (|48p into a one-dimensional effective equation along the axial direction. 
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First, integrating (j28b|) over i?(x) for x G dO, we obtain the following no-flux 
boundary condition: 

[VxP - f (x)p] ■ n = on dO. (49) 

Analogously to the point-particles case, we use the narrow-domain variables (fS]) 
and define p(x, t) = ep(x, t). With this rescaling, Eqs. (|48(1 and (|49|) become 



|?=/ip on i = ±i, (50b) 
ax 2 

|| = e/2P on = (50c) 

where fi = fi{x, ey). There is no need to expand the integral terms X^*-* in terms of 
the narrow-domain variables since these are written in terms of the inner variables 
(|35p , which coincide with the narrow-domain variable for expressions independent 
of the second particle's coordinates. 

Expanding p in powers of e and solving (|5Up gives, at leading order, that p is 
independent of y. As before, we introduce the effective one-dimensional densities 
as pe*^ = /^02P^''''^3^- Thus we have that pi°^ = hpi'^K At the next order, 

p(i)(x,t) = f2{x,0)p^°\x,t)y + pi^\x,t)/h. (51) 

For clarity of notation, in the remaining of this section we write /,;(x,0) = fi. 
Integrating the second order of (|50ap over the channel's cross section and using 
(|50cp . gives 

where we have used that f^[^^^l^^My = [see Eqs. (gS]) and (gTa])]. Note that this 
equation coincides with the effective equation for point particles ([7]). It is at the 
next order that the finite-size effects appear. 

Repeating the same procedure of integrating with respect to y the 0{e^) of 
(|50p and using (|50cp to eliminate p^^-* yields the following solvability condition 



dt 



Using (|46p and (|47p . the cross-section integral of I*-^-* is 

h/2 ^1 ^ \ y / ^g^^ 
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where Mi(h) = J^j^^^ iii{h, y)dy reads 
Mi{h) =TTh-^+0{l-h) 



- (2 + /i^)a/1 - ft2 _ 2h arccos(/i) 
o 



where 0{x) is the Heaviside step-function, and M[Q\{h) = /^^j ^W- 

Although Qi and Q2 are only solved numerically, using information from their 
respective problems (|4ip and (|42|) one can deduce analytical expressions for their 
integrals M\ Qj], namely t hat M[Qi] = -Mi(/i)/(2\/2) and M[Q2] = (see Ap- 
pendix C.3 in lBruna|[2012h . Using this, we find that 



(55) 



h/2 



Because q is independent of the inner variables (i, 3/1, j/2), we can write q{xi,t) = 
qe{x, t)/h. Moreover, the normalisation condition on P gives that qe{x, t) = pe{x, t)+ 
Oie). Therefore the right-hand side of ([56]) becomes ^(pe^)/ft^. 
Combining dSH), ^ and ^ yields 



which coincides with the effective equation (|10|) for a two-dimensional narrow- 
channel after writing a/j = Mi{h)/h^; see (|lip . 



5.3 Outline of steps for other geometries 

In this section we indicate the key steps to derive the effective continuum Fokker- 
Planck equation (|10|) for a general geometry, and in particular for the three- 
dimensional cases (NC3) and (PP) presented in Section [2.31 First, we note below 
relevant definitions that change with the problem dimension d and the number of 
confined dimensions k (recall that de = d — k): 

(1) Identify the number of confined and effective dimensions: the original position 
vector is split into two components, x = (xe,Xc) (see Table [ij. For example, for 
(NC3) Xe = X and Xc = {y,z), while for (PP) Xe = {x,y) and Xc = z. 



Dimension Variables Domain 

Original problem d x i7 

Confined space k Xc f^c 

Effective problem de Xe fie 

Table 1 Dimension and variables in each of the spaces and subspaces. 



(2) Determine the confinement parameter (s): Next we must choose an scaling for 
the confined dimensions relative to the unconfined ones. In all cases considered 
here we made that simple by saying that all confined dimensions are of length H 
relative to the unconfined ones, but there could be of different lengths too, such as 
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the narrow channel of rectangular cross section mentioned briefly in Subsec. 13.31 
Suppose that the confinement dimensions are H = (Hi, . . . , Hk) = 0{e). Then the 
vector of confinement parameters is given hy h = {hi, . . . , ftj.), with = H^/e. 
Note that we are assuming that a confined dimension is always of order e (the 
particle's diameter). However, this could also be generalised and introduce an 
intermediate scaling (between e and one) . 

(3) Narrow-domain variables transformation: the generalisation of the change of 
variables ([5]) is 

Xe — Xe, Xc — eXc- 

We write x = {xe,Xc)- The one-particle and two-particle densities in the rescaled 
domain are respectively defined as 

p{xi,t) = e''p(xi,t), P(xi,X2,f) = e^''P(xi,X2,t). 

The original domain O is mapped into the rescaled domain cj. 

(4) Apply the effective domain transformation: The 72 transformation to reduce the 
original d-dimensional problem to a de-dimensional problem (cf. Subsec. 15. 2p re- 
quires the following rescaled densities 

Pe{xi,t) = ylp(xi,t), Pe{xi,X2,t) = yl^P(xi,X2,t), 

where 

^ = n■=lh^■ (58) 

This factor is, in fact, equal to the volume of the rescaled domain oj (also = |a;c|). 
Recall it is introduced so that pe and Pe are defined as densities. 

(5) Evaluate the collision integral I: The evaluation of the contribution of the two- 
particle interaction I reduces to computing one coefficient like M\(h) in (|57|) for 
each of the unconfined dimensions. Suppose that all the unconfined dimensions 
are symmetric. In general it is equal to 

Mi{h)= [ pi{h,xc)dxc, (59) 

where Xc are the confined coordinates of the first particle (it was simply y in 
the (NC2) [cf. (|55|l ]. The function pi(h,Xc) is the integral of Sy^ over the contact 
surface between two particles when the first one has coordinates (xe,Xc). Without 
loss of generality we can set Xe = Oe. In the rescaled problem, this surface is a 
d— dimensional unit sphere centred at {Oe,Xc), and (possibly) intersected with the 
confinement walls dujc'- 

tii{h,xc)= I x^dS, (60) 

where B{xc) is the unit ball, dS is the surface differential and x is the unconfined 
dimension of this surface. Once AIi is computed, we use it for the nonlinear coeffi- 
cient of the effective Fokker-Planck equation. The generalisation of the coefficient 
Ofi in (|10p is given by 

ah = ^Mi{h). (61) 
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6 Discussion 



In this paper, we have presented a derivation of a continuum model for the diffusion 
of finite-size particles in a confined domain whose dimensions are comparable to 
the particle dimensions. We have given the model explicitly for three confined 
geometries, namely a two- and three-dimensional narrow channel of square cross- 
section and Hele-Shaw cell (two close parallel plates), but also indicated how the 
model derivation for more general geometries can be done. The resulting continuum 
model is a nonlinear drift-diffusion equation for the one-particle probability density 
or for the population volume concentration, with the nonlinear term depending on 
the excluded- volume created by the particles as well as the confinement parameter 
(in the cases we have studied, the confinement parameter is the channel width or 
the separation between plates). This equation is defined in the effective domain, a 
lower-dimensional domain for the unconfined dimensions only, exploiting the fact 
that equilibration in the confined dimensions is relatively fast and most dynamics 
occur along the unconfined ones. For example, in the case of a two- or three- 
dimensional channel the resulting continuum model is a one-dimensional equation 
along the axial direction. 

The derivation of the final continuum model involv ed two key steps. First, 
as in our previous works (jBruna and Chapman|[2012allbl) . we used the method 
of matched asymptotic expansions to reduce the high-dimensional Fokker-Planck 
equation associated with the individual-based description of the system to a low- 
dimensional Fokker-Planck equation for the one-particle density function. Second, 
we exploited the confined geometry of the domain to perform a further reduction 
and integrate out the confined dimensions of the problem. While in our previous 
works the particle— particle-wall interactions were a higher-order correction, and 
thus were neglected, in this work they were taken into account. In other words, 
while in an unconfined domain such three-body interactions simply create a bound- 
ary layer near the domain walls, in the situations considered in the present work 
this boundary layer extends across the cross section and must be solved accurately. 

The model has two interesting features. First, we found that, for a given volume 
fraction, there exists an optimal ratio h between confinement and particle size such 
that the excluded-volume effects are maximised. This means, for example, that one 
can design a lab-on-a-chip device to achieve a maximal collective diffusion coeffi- 
cient. Second, the model is capable of describing the whole range of confinement 
levels and interpolating between confinement extremes, that is, between extreme 
confinement {h = 0) when particles cannot pass each other to unconfined diffusion 
{h = oo). We have examined the limiting models corresponding to each of the 
three geometries we have studied, and found that our model agrees with them. 
For example, in the narrow-ch annel cases , the extreme confinement corresponds 
to a single-file diffusion model (|Rostlll984l) while the other limit is an unconfined 
two or three-dimensional diffusion model ( Bruna and Chapman|[2012bl ). This is a 
useful analytical tool to predict the error that is being committed by using the 
limiting models, that is, either ignoring the fact that particles can (quite) pass 
each other and using the purely one-dimensional single-file equation, or neglecting 
confinement conditions and boundary layers. 

In order to assess the validity of our model predictions, we have compared the 
numerical results of our continuum model for the two-dimensional narrow channel 
to the results from the corresponding individual-based model as well as the limiting 
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continuum models of point particles, single-file diffusion and unconfined diffusion. 
We have observed excellent agreement between the narrow channel model and 
the stochastic simulations of the discrete model, and confirmed the interpolating 
properties of the model between confinement extremes. Finally, we have discussed 
a case study involving the diffusion under a ratchet potential, used for example to 
describe the transport of molecular motors through microtubules. We have exam- 
ined the effects that interactions between particles and confinement conditions can 
have in the analysis of the problem. For example, we have found that an increase 
on excluded-volume effects causes the particle transport to be less sensitive to the 
tilting of the ratchet potential. 

There are several directions along which one can extend this model. Firstly, it 
would be interesting to allow mul tiple types of particles (labelle d blues and reds, 
say) by combining the analysis in iBruna and ChapmanI ( 2012a|) with the present 
work. A challenging aspect of this problem is that the resulting model of two 
species in a narrow channel must account for the fact that as the channel becomes 
single file the order of particles is fixed, and since the two types of particles are 
distinguishable the order of particles matters. In other words, we expect a qualita- 
tive difference in the model as h cross the value of one: for h > 1, even if all the red 
particles were initially to the left of the blue particles, we expect the two popula- 
tions to mix together (that is, the jump in the initial density profiles will spread); 
in contrast, for h < 1, if the two populations are segregated initially, they will stay 
like that for all times, and we expect to see this refelected in the population-level 
densities. In connection with this issue is the fac t that the self-diffusion of a par - 
ticle is not defined in one-dimensional systems ( Ackerson and Fleishmanlll982l) . 
Another extension is to allow channels of varying cross section, that is, to extend 
the Fick-Jacobs model ^ to the case of finite-size particles. A simplification of 
this problem (igno ring the particle-particle interactions) has been considered by 
iRiefier et all ( 2010l) . A related problem would be to match a narrow channel with 
two bulk domains in each end; such a geometry could have important applications 
in the area of ion channels. 
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